home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / modelers / geomview / source.lha / Geomview / src / lib / geometry / hpointn / hpointnb.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-09-03  |  7.8 KB  |  350 lines
  1. /* Copyright (c) 1993 The Geometry Center; University of Minnesota
  2.    1300 South Second Street;  Minneapolis, MN  55454, USA;
  3.    
  4. This file is part of geomview/OOGL. geomview/OOGL is free software;
  5. you can redistribute it and/or modify it only under the terms given in
  6. the file COPYING, which you should have received along with this file.
  7. This and other related software may be obtained via anonymous ftp from
  8. geom.umn.edu; email: software@geom.umn.edu. */
  9.  
  10. /* Authors: Charlie Gunn, Pat Hanrahan, Stuart Levy, Tamara Munzner, Mark Phillips */
  11.  
  12.  
  13. #include <ooglutil.h>
  14. #include "hpointn.h"
  15. #include <math.h>
  16.  
  17.  
  18. HPointN *
  19. HPtNCopy(const HPointN *pt1, HPointN *pt2)
  20. {
  21.     if(pt2 == NULL) {
  22.         pt2 = HPtNCreate(pt1->dim, pt1->v);
  23.     } else {
  24.         if(pt2->dim != pt1->dim) {
  25.             pt2->v = OOGLRenewNE(HPtNCoord, pt2->v, pt1->dim, "renew HPointN");
  26.             pt2->dim = pt1->dim;
  27.         }
  28.         memcpy(pt2->v, pt1->v, pt1->dim*sizeof(HPtNCoord));
  29.     }
  30. /*    pt2->space = pt1->space; */
  31.     return(pt2);
  32. }
  33.  
  34. HPointN *
  35. HPtNPad(HPointN *pt1, short dim2, HPointN *pt2)
  36. {
  37.     short dim1 = pt1->dim;
  38.     int i;
  39.  
  40.     if( dim2 < 1 ) {
  41.         return(NULL);
  42.     }
  43.  
  44.     if(pt1 == NULL) {
  45.         pt2 = HPtNCreate(dim2,NULL);
  46.         return(pt2);
  47.     }
  48.  
  49.  
  50.     if(pt1 != pt2) {
  51.         if(pt2 == NULL)
  52.             pt2 = HPtNCreate(dim2,NULL);
  53.         else if (pt2->dim != dim2) {
  54.                 pt2->v = OOGLRenewNE(HPtNCoord, pt2->v, dim2, "renew HPointN");
  55.             pt2->dim = dim2;
  56.         }
  57.         if (dim1 <= dim2) {
  58.             bzero(pt2->v,dim2*sizeof(HPtNCoord));
  59.             for(i=0; i<dim1-1; i++)
  60.                 pt2->v[i] = pt1->v[i];
  61.         } else if ( dim1 > dim2 ) {
  62.             for(i=0; i<dim2-1; i++)
  63.                 pt2->v[i] = pt1->v[i];
  64.         }
  65.         pt2->v[dim2-1] = pt1->v[dim1-1];
  66.     } else {
  67.         HPointN *newpt = HPtNCreate(dim2,NULL);
  68.         if (dim1 <= dim2) {
  69.             for(i=0; i<dim1-1; i++)
  70.                 newpt->v[i] = pt1->v[i];
  71.         } else if ( dim1 > dim2 ) {
  72.             for(i=0; i<dim2-1; i++)
  73.                 newpt->v[i] = pt1->v[i];
  74.         }
  75.         newpt->v[dim2-1] = pt1->v[dim1-1];
  76.         HPtNCopy(newpt,pt2);
  77.     }
  78.     return(pt2);
  79. }
  80.  
  81. HPointN *
  82. HPtNAdd(const HPointN *pt1, const HPointN *pt2, HPointN *sum)
  83. {
  84.     short dim1 = pt1->dim, dim2 = pt2->dim;
  85.     int i;
  86.     
  87.   if(dim1 == dim2) {
  88.       if(sum == NULL)
  89.         sum = HPtNCreate(dim1,NULL);
  90.     else if(sum->dim < dim1) {
  91.            sum->v = OOGLRenewNE(HPtNCoord, sum->v, dim1, "renew HPointN");
  92.         sum->dim = dim1;
  93.     }
  94.     for( i=0; i<dim1-1; i++)
  95.         sum->v[i] = pt1->v[i] + pt2->v[i];
  96.     sum->v[dim1-1] = pt1->v[dim1-1];
  97.     sum->dim = dim1;
  98.   } else if (dim1 > dim2) {
  99.       if(sum == NULL)
  100.         sum = HPtNCreate(dim1,NULL);
  101.     else if(sum->dim < dim1) {
  102.            sum->v = OOGLRenewNE(HPtNCoord, sum->v, dim1, "renew HPointN");
  103.         sum->dim = dim1;
  104.     }
  105.     for( i=0; i<dim2-1; i++)
  106.         sum->v[i] = pt1->v[i] + pt2->v[i];
  107.     for( i=dim2-1; i<dim1-1; i++)
  108.         sum->v[i] = pt1->v[i];
  109.     sum->v[dim1-1] = pt1->v[dim1-1];
  110.     sum->dim = dim1;
  111.   } else if (dim1 < dim2) {
  112.       if(sum == NULL)
  113.         sum = HPtNCreate(dim2,NULL);
  114.     else if(sum->dim < dim2) {
  115.            sum->v = OOGLRenewNE(HPtNCoord, sum->v, dim2, "renew HPointN");
  116.         sum->dim = dim2;
  117.     }
  118.     for( i=0; i<dim1-1; i++)
  119.         sum->v[i] = pt1->v[i] + pt2->v[i];
  120.     for( i=dim1-1; i<dim2-1; i++)
  121.         sum->v[i] = pt2->v[i];
  122.     sum->v[dim2-1] = pt1->v[dim1-1];
  123.     sum->dim = dim2;
  124.   }
  125.   return(sum);
  126.  
  127. }
  128.     
  129.     /* Space */
  130. int 
  131. HPtNSpace( const HPointN *pt );
  132. HPointN *
  133. HPtNSetSpace( HPointN *pt, int space );
  134.  
  135.     /* Linear combination */
  136. HPointN *
  137. HPtNComb(HPtNCoord u, const HPointN *pu, HPtNCoord v, const HPointN *pv, HPointN *sum)
  138. {
  139.     short dimu = pu->dim, dimv = pv->dim;
  140.     int i;
  141.  
  142.   if(dimu == dimv) {
  143.       if(sum == NULL)
  144.         sum = HPtNCreate(dimu,NULL);
  145.     else if(sum->dim < dimu) {
  146.            sum->v = OOGLRenewNE(HPtNCoord, sum->v, dimu, "renew HPointN");
  147.         sum->dim = dimu;
  148.     }
  149.     for( i=0; i<dimu-1; i++)
  150.         sum->v[i] = u * pu->v[i] + v * pv->v[i];
  151.     sum->v[dimu-1] = pu->v[dimu-1];
  152.     sum->dim = dimu;
  153.   } else if( dimu > dimv ) {
  154.       if(sum == NULL)
  155.         sum = HPtNCreate(dimu,NULL);
  156.     else if(sum->dim < dimu) {
  157.            sum->v = OOGLRenewNE(HPtNCoord, sum->v, dimu, "renew HPointN");
  158.         sum->dim = dimu;
  159.     }
  160.     for( i=0; i<dimv-1; i++)
  161.         sum->v[i] = u * pu->v[i] + v * pv->v[i];
  162.     for( i=dimv-1; i<dimu-1; i++)
  163.         sum->v[i] = u * pu->v[i];
  164.     sum->v[dimu-1] = pu->v[dimu-1];
  165.     sum->dim = dimu;
  166.   } else if( dimu < dimv ) {
  167.       if(sum == NULL)
  168.         sum = HPtNCreate(dimv,NULL);
  169.     else if(sum->dim < dimv) {
  170.            sum->v = OOGLRenewNE(HPtNCoord, sum->v, dimv, "renew HPointN");
  171.         sum->dim = dimv;
  172.     }
  173.     for( i=0; i<dimu-1; i++)
  174.         sum->v[i] = u * pu->v[i] + v * pv->v[i];
  175.     for( i=dimu-1; i<dimv-1; i++)
  176.         sum->v[i] = v * pv->v[i];
  177.     sum->v[dimv-1] = pu->v[dimu-1];
  178.     sum->dim = dimv;
  179.   }
  180.   return(sum);
  181. }
  182.  
  183.     /* Reduce to unit vector */
  184. HPointN *
  185. HPtNUnit(const HPointN *from, HPointN *to)
  186. {
  187.     register HPtNCoord len = 0;
  188.     short dim = from->dim;
  189.     int i;
  190.  
  191.     for( i=0; i<dim-1; i++)
  192.         len += from->v[i] * from->v[i];
  193.  
  194.     if( (len = sqrt( len )) != 0) {
  195.           if(to == NULL)
  196.             to = HPtNCreate(dim,NULL);
  197.         else if(to->dim != dim) {
  198.                to->v = OOGLRenewNE(HPtNCoord, to->v, dim, "renew HPointN");
  199.             to->dim = dim;
  200.         }
  201.         for( i=0; i<dim-1; i++)
  202.             to->v[i] = from->v[i] /len;
  203.         to->v[dim-1] = from->v[dim-1];
  204.     } else
  205.         HPtNCopy(from,to);
  206.  
  207.     return(to);
  208. }
  209.  
  210.     /* Dehomogenize */
  211. HPtNCoord
  212. HPtNDehomogenize(const HPointN *from, HPointN *to)
  213. {
  214.     short dim = from->dim;
  215.     register HPtNCoord inv = 1.0 / from->v[dim-1];
  216.     int i;
  217.  
  218.     
  219.     if( from->v[dim-1] == 1.0 || from->v[dim-1] == 0.0) {
  220.         if( from != to) HPtNCopy(from,to);
  221.         return;
  222.     }
  223.  
  224.  
  225.       if(to == NULL)
  226.         to = HPtNCreate(dim,NULL);
  227.     else if(to->dim != dim) {
  228.           to->v = OOGLRenewNE(HPtNCoord, to->v, dim, "renew HPointN");
  229.         to->dim = dim;
  230.     }
  231.  
  232.     for( i=0; i<dim-1; i++)
  233.         to->v[i] = from->v[i] * inv;
  234.     to->v[dim-1] = 1.0;
  235.     return(1./inv);
  236. }
  237.     
  238.  
  239.     /* Apply a TransformN to an HPointN */
  240. HPointN *
  241. HPtNTransform( const TransformN *T, const HPointN *from, HPointN *to )
  242. {
  243.     short idim = T->idim, odim = T->odim, dim = from->dim;
  244.     int i, j;
  245.     HPtNCoord *v = (HPtNCoord *)malloc(idim*sizeof(HPtNCoord));
  246.  
  247.     for( i=0; i<idim; i++)
  248.         v[i] = from->v[i];
  249.  
  250.       if(to == NULL)
  251.         to = HPtNCreate(odim,NULL);
  252.     else if(to->dim != odim) {
  253.           to->v = OOGLRenewNE(HPtNCoord, to->v, odim, "renew HPointN");
  254.         to->dim = odim;
  255.     }
  256.  
  257.   if (idim == dim) {
  258.     for( i=0; i<odim; i++) {
  259.         to->v[i] = 0;
  260.         for( j=0; j<idim; j++)
  261.             to->v[i] += v[j] * T->a[j*odim+i];
  262.     }
  263.   } else if (idim > dim) { /* padding with zeroes */
  264.     for( i=0; i<odim; i++) {
  265.         to->v[i] = 0;
  266.         for( j=0; j<dim; j++)
  267.             to->v[i] += v[j] * T->a[j*odim+i];
  268.     }
  269.   } else if (idim < dim) {
  270.     for( i=0; i<odim; i++) {
  271.         to->v[i] = 0;
  272.         for( j=0; j<idim; j++)
  273.             to->v[i] += v[j] * T->a[j*odim+i];
  274.         if( i>idim-1)
  275.             to->v[i] += v[i];
  276.     }
  277.   }
  278.  
  279.   free(v);
  280.   return(to);
  281. }
  282.  
  283.     /* Dot product of two vectors */
  284. HPtNCoord 
  285. HPtNDot( const HPointN *p1, const HPointN *p2)
  286. {
  287.     HPtNCoord result;
  288.     register int i;
  289.     short dim = p1->dim;
  290.  
  291.     if (p2->dim < dim) {
  292.         dim = p2->dim;
  293.     }
  294.  
  295.     result = 0;
  296.     for( i = 0; i< dim-1; i++)
  297.         result += p1->v[i] * p2->v[i];
  298.  
  299.     return(result);
  300. }
  301.  
  302.     /* Return index'th component of p . T */
  303.     /* If index is out of range (e.g. -1), return N-1'th component
  304.      * (i.e. homogeneous divisor).
  305.      */
  306. HPtNCoord *
  307. HPtNTransformComponents( const HPointN *p, const TransformN *T, int ncomponents, int *indices, HPtNCoord *results )
  308. {
  309.     short idim = T->idim, odim = T->odim, dim = p->dim;
  310.     int i, j, k;
  311.     HPtNCoord *v = p->v;
  312.  
  313.     if( ncomponents < 0 || ncomponents > odim ) {
  314.         return NULL;
  315.     }
  316.  
  317.     if ( idim == dim ) {
  318.          for ( k = 0; k < ncomponents; k++ ) {
  319.             i = indices[k];
  320.             if( i < 0 || i > odim )
  321.                 i = odim-1;
  322.             results[k] = 0;
  323.             for( j=0; j<idim; j++)
  324.                 results[k] += v[j] * T->a[j*odim+i];
  325.         }
  326.     } else if ( idim < dim ) {
  327.          for ( k = 0; k < ncomponents; k++ ) {
  328.             i = indices[k];
  329.             if( i < 0 || i > odim )
  330.                 i = odim-1;
  331.             results[k] = v[dim-1] * T->a[(idim-1)*odim + i];
  332.             for( j=0; j<idim-1; j++)
  333.                 results[k] += v[j] * T->a[j*odim+i];
  334.             if( i>idim-1)
  335.                 results[k] += v[i];
  336.         }
  337.     } else { /* idim > dim */
  338.          for ( k = 0; k < ncomponents; k++ ) {
  339.             i = indices[k];
  340.             if( i < 0 || i > odim )
  341.                 i = odim-1;
  342.             results[k] = v[dim-1] * T->a[(idim-1)*odim + i];
  343.             for( j=0; j<dim-1; j++)
  344.                 results[k] += v[j] * T->a[j*odim+i];
  345.         }
  346.     }
  347.  
  348.   return(results);
  349. }
  350.